%ֻ������ֵ��Ϊ������������
%clc;
%clear all;   
%iII=1;
function thresholding(num)

iI=num;
t=5;
file='test';
file1=strcat(file,num2str(iI));
p1=strcat(file1,'.tif');

file2=strcat(file,num2str(iI+1));
p2=strcat(file2,'.tif');

t_k5=cputime;
I0=imread(p1);
[Hh,Ww]=size(I0);
%�Աȶ���ǿ
I=ctrenhance(p1);
%I=double(I);
x21=mat2gray(I);%ת��Ϊ�Ҷ�ͼ�� 0-1
[K11,K21]=size(x21);

%%%%%%%%%%%%%%%%��ֵ�˲��õ�y1%%%%%%%%%%%%%%%%%%
m=3;
I11=medfilt2(x21,[m,m]); 
% figure,imshow(I11)
%��I��������λ��para1����һЩ��������morphspotgrid1��ע��
%para2Ϊ����ˮƽ������λ�õľ���
%para3Ϊ���ش�ֱ������λ�õľ���
%hnum=para1(5);%ˮƽ������
%vnum=para1(6);%��ֱ������
%[para1,para2,para3]=morphspotgrid1(I);
[para1,para2,para3,para4,para5]=morphspotgrid(I11);
para2=round(para2);
if para2(1)==0
    para2(1)=1;
end
para3=round(para3);
if para3(1)==0
    para3(1)=1;
end


hnum=size(para2,2);%ˮƽ������
vnum=size(para3,2);%��ֱ������

% para2
I=uint8(I11.*256);
[h,w]=size(I);
h1=size(para4,2);
v1=size(para5,2);
% hnum=hnum-3;
% h1=h1-5;
% para4=para4(1:h1);
% hnum=hnum-1;
% h1=h1-2;
% para4(1:11)=para41(1:11);
% para4(12:h1)=para41(14:h1+2);
% hnum=hnum+1;
% h1=h1+2;
% para4(1:h1-2)=para41(1:h1-2);
% para4(h1-1)=para4(h1-2)+10;
% para4(h1)=para4(h1-1)+16;
% para2(hnum-1)=para2(hnum-2)+12;
% para2(hnum)=398;
t_k51=cputime;
% t=4;
 
for j=1:h
    if j<para4(2)-t
        I(j,:)=0;
    end
    for i=3:2:h1-1
       if j>para4(i)+t&&j<para4(i+1)-t
         I(j,:)=0;
       end
    end
    if j>para4(h1)+t
        I(j,:)=0;
    end
end


for j=1:w
    if j<para5(2)-t
       I(:,j)=0;
    end
    for i=1:2:v1-1
      if j>para5(i)+t&&j<para5(i+1)-t
        I(:,j)=0;
      end
    end
    if j>para5(v1)+t
        I(:,j)=0;
    end
end

% figure,imhist(I),title('����оƬͼ��ĻҶ�ֱ��ͼ');
% [X1,map]=imhist(I);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
III12=zeros(h,w);
for i1=1:hnum-1
    for j1=1:vnum-1        
%         for i2=1:para2(i1+1)-para2(i1)+1  
%             for j2=1:para3(j1+1)-para3(j1)+1
%                 I1(i2,j2)=I(i2+para2(i1)-1,j2+para3(j1)-1); 
        for i2=1:para4(2*i1+1)-para4(2*i1)+2*t-1
            for j2=1:para5(2*j1+1)-para5(2*j1)+2*t-1
                I1(i2,j2)=I(i2+para4(2*i1)-t,j2+para5(2*j1)-t);
            end
        end
        [H,W]=size(I1);
[X1,map]=imhist(I1);
ii=0;
num=H*W;
II=zeros(num,1);
%����������ֵ����������II
for i=1:H
    for j=1:W
        ii=ii+1;
        II(ii,1)=I1(i,j);
    end
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%��yu1Ϊ��ֵ����ÿ��С����Ļ����ֳ�ǰ���ͱ���
s=0;
for i=1:max(map)+1
    s=s+X1(i)*(i-1);
end
yu1=s/sum(X1);
yu2=median(median(I1));
iii=0;
for i=1:H
    for j=1:W
        iii=iii+1;
        if II(iii)>yu1
%             III12(i+para2(i1)-1,j+para3(j1)-1)=1;
            III12(i+para4(2*i1)-t,j+para5(2*j1)-t)=1;
        end
    end
end
 clear I1
    end
end

% figure,imshow(III12),title('k_cy3');
% imwrite(III12,'k_cy1.tif','tif');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%�󲹳�����ԭͼ�Ͻ��ָ��Ե����
%kmeans����󻮷ֱ�Ե
k_edge3=edge(III12,'canny');

PP3=imread(p1);
[h1,w1]=size(PP3);
figure,imshow(PP3)
% set(gca,'pos',[0 0 1 1],'DataAspectRatioMode','auto')
% set(gca,'border','tight','initialmagnification','fit');
% set(gca,'looseInset',[0 0 0 0])
% set(gca,'LooseInset',get(gca,'TightInset'))

for i=1:h1
  for j=1:w1
         if k_edge3(i,j)==1
            hold on
         plot(j,i,'r');
         end
  end
end
%title('edge of p1 ');
time_k5=cputime-t_k5
time_k51=cputime-t_k51
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%�������յõ��ķָ�ͼ������������ݵĸ�������
%������������ڵĸ�������
PP1=imread(p1);
PP1=double(PP1);
PP1=PP1./256;        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%��ǰ���������صĻҶ�ֵ��CH1_RAW�ͱ����������صĻҶ�ֵ��CH1_BKD_RAW
%��ǰ���������ҶȾ�ֵ����CH1_MEAN-CH1_BKD_MEAN=CH1_MINUS,ǰ���ͱ����ĻҶȾ�ֵ��
%ͬʱ�����������ڸ����������ǰ�����ظ���CH1_AREA�ͱ������ظ���CH1_BKD_AREA

%CH1_MEAN_MINUS=zeros(hnum-1,vnum-1);%ǰ���ͱ����ĻҶȾ�ֵ��
CH1_RAW=zeros(hnum-1,vnum-1);%ǰ���������صĻҶ�ֵ��
CH1_BKD_RAW=zeros(hnum-1,vnum-1);%�����������صĻҶ�ֵ��
CH1_AREA=zeros(hnum-1,vnum-1);%��¼(i,j)������ ǰ�����ص����
CH1_BKD_AREA=zeros(hnum-1,vnum-1);%��¼(i,j)������ �������ص����
CH1_MEAN=zeros(hnum-1,vnum-1);%ǰ���ĻҶȾ�ֵ
% CH2_MEAN=zeros(hnum-1,vnum-1);
CH1_BKD_MEAN=zeros(hnum-1,vnum-1);%�����ĻҶȾ�ֵ

for i=1:hnum-1
    for j=1:vnum-1
       for ii=para2(i):para2(i+1)
         for jj=para3(j):para3(j+1)
%         for ii=para4(2*i)-t+1:para4(2*i+1)+t-1
%            for jj=para5(2*j)-t+1:para5(2*j+1)+t-1
              %if ii<=h&&jj<=w
            if (III12(ii,jj)==1)
               CH1_RAW(i,j)=CH1_RAW(i,j)+PP1(ii,jj);%ǰ���ĻҶ�ֵ��
               CH1_AREA(i,j)=CH1_AREA(i,j)+1;%ǰ�������ظ���
            else
            
            CH1_BKD_RAW(i,j)=CH1_BKD_RAW(i,j)+PP1(ii,jj);%�����ĻҶ�ֵ��
            CH1_BKD_AREA(i,j)=CH1_BKD_AREA(i,j)+1;%���������ظ���
            end
              %end
          end
        end
if CH1_BKD_AREA(i,j)==0    
   CH1_BKD_MEAN(i,j)=0;%�󱳾��ĻҶȾ�ֵ
else
   CH1_BKD_MEAN(i,j)=CH1_BKD_RAW(i,j)/CH1_BKD_AREA(i,j);%�󱳾��ĻҶȾ�ֵ  
end
% CH1_BKD_MEAN(i,j)=CH1_BKD_RAW(i,j)/CH1_BKD_AREA(i,j);%�󱳾��ĻҶȾ�ֵ
if CH1_AREA(i,j)==0    
   CH1_RAW(i,j)=1; 
   CH1_MEAN(i,j)=0;%��ǰ���ĻҶȾ�ֵ
else
   CH1_MEAN(i,j)=CH1_RAW(i,j)/CH1_AREA(i,j);%��ǰ���ĻҶȾ�ֵ
  
end
 CH1_MEAN(i,j)=CH1_MEAN(i,j)-CH1_BKD_MEAN(i,j);
    end
end
% xlswrite('k_CH1_AREA.xls',CH1_AREA)
% xlswrite('k_CH1_BKD_AREA.xls',CH1_BKD_AREA)
% xlswrite('k_CH1_MEAN.xls',CH1_MEAN)
% xlswrite('k_CH2_MEAN.xls',CH2_MEAN)
% xlswrite('k_CH1_BKD_MEAN.xls',CH1_BKD_MEAN)
% xlswrite('k_CH1_RAW.xls',CH1_RAW)
% xlswrite('k_CH1_BKD_RAW.xls',CH1_BKD_RAW)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%��ɹ���
flag=0;
medi1=median(median(CH1_AREA));
for i=1:hnum-1
    for j=1:vnum-1
        if CH1_AREA(i,j)>medi1*0.4&&CH1_AREA(i,j)<medi1*1.5
            flag=flag+1;
        end
    end
end
success1=flag/(i*j)
% success1


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t_k5=cputime;
%ֻ������ֵ��Ϊ������������
I02=imread(p2);
[Hh,Ww]=size(I02);
%�Աȶ���ǿ
I2=ctrenhance(p2);
%I2=double(I2);
x22=mat2gray(I2);%ת��Ϊ�Ҷ�ͼ�� 0-1
[K12,K22]=size(x22);

%%%%%%%%%%%%%%%%��ֵ�˲��õ�y1%%%%%%%%%%%%%%%%%%
m=3;
I21=medfilt2(x22,[m,m]); 
%��I2��������λ��para21����һЩ��������morphspotgrid1��ע��
%para22Ϊ����ˮƽ������λ�õľ���
%para23Ϊ���ش�ֱ������λ�õľ���
%hnum2=para21(5);%ˮƽ������
%vnum2=para21(6);%��ֱ������
%[para21,para22,para23]=morphspotgrid1(I);

[para1,para2,para3,para4,para5]=morphspotgrid(I21);
para2=round(para2);
if para2(1)==0
    para2(1)=1;
end
para3=round(para3);
if para3(1)==0
    para3(1)=1;
end


hnum=size(para2,2);%ˮƽ������
vnum=size(para3,2);%��ֱ������
%para4
%para5

I2=uint8(I21.*256);
[h,w]=size(I2);
h2=size(para4,2);
v2=size(para5,2);
% para2
% para4

% hnum=hnum-3;
% h2=h2-5;
% para4=para4(1:h2);
% hnum=hnum-1;
% h2=h2-2;
% para4(1:h2)=para41(1:h2);
% para4(h2)=para4(h2)+10;
t_k51=cputime;
% t=4;

for j=1:h
    if j<para4(2)-t
        I2(j,:)=0;
    end
    for i=3:2:h2-1
       if j>para4(i)+t&&j<para4(i+1)-t
         I2(j,:)=0;
       end
    end
    if j>para4(h2)+t
        I2(j,:)=0;
    end
end


for j=1:w
    if j<para5(2)-t
        I2(:,j)=0;
    end
    for i=1:2:v2-1
      if j>para5(i)+t&&j<para5(i+1)-t
        I2(:,j)=0;
      end
    end
    if j>para5(v2)+t
        I2(:,j)=0;
    end
end

% figure,imhist(I2),title('imhist of p2');
% [X2,map2]=imhist(I2);

III22=zeros(h,w);
for i1=1:hnum-1
    for j1=1:vnum-1        
%         for i2=1:para2(i1+1)-para2(i1)+1  
%             for j2=1:para3(j1+1)-para3(j1)+1
%                 I12(i2,j2)=I(i2+para2(i1)-1,j2+para3(j1)-1); 
        for i2=1:para4(2*i1+1)-para4(2*i1)+2*t-1
            for j2=1:para5(2*j1+1)-para5(2*j1)+2*t-1
                I12(i2,j2)=I2(i2+para4(2*i1)-t,j2+para5(2*j1)-t);
            end
        end
        [H,W]=size(I12);
% [H,W]=size(I2);
[X2,map2]=imhist(I12);
ii=0;
num=H*W;
II2=zeros(num,1);
%����������ֵ����������II
for i=1:H
    for j=1:W
        ii=ii+1;
        II2(ii,1)=I12(i,j);
    end
end
clear I12

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%��yu2Ϊ��ֵ����ÿ��С����Ļ����ֳ�ǰ���ͱ���
s=0;
for i=1:max(map2)+1
    s=s+X2(i)*(i-1);
end
yu2=s/sum(X2);
iii=0;
for i=1:H
    for j=1:W
        iii=iii+1;
        if II2(iii)>yu2
%             III22(i+para2(i1)-1,j+para3(j1)-1)=1;
            III22(i+para4(2*i1)-t,j+para5(2*j1)-t)=1;
        end
    end
end
    end
end


% figure,imshow(III22),title('k_cy5');
% imwrite(III22,'k_cy2.tif','tif');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%�󲹳�����ԭͼ�Ͻ��ָ��Ե����
%kmeans����󻮷ֱ�Ե
k_edge5=edge(III22,'canny');

PP5=imread(p2);
[h1,w1]=size(PP5);
figure,imshow(PP5,[]);
for i=1:h1
  for j=1:w1
         if k_edge5(i,j)==1
            hold on
         plot(j,i,'r');
         end
  end
end
% title('edge of p2 ');
time_k5=cputime-t_k5
time_k51=cputime-t_k51
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%�������յõ��ķָ�ͼ������������ݵĸ�������
%������������ڵĸ�������
        
PP2=imread(p2);
PP2=double(PP2);
PP2=PP2./256;  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%��ǰ���������صĻҶ�ֵ��CH2_RAW�ͱ����������صĻҶ�ֵ��CH2_BKD_RAW
%��ǰ���������ҶȾ�ֵ����CH2_MEAN-CH2_BKD_MEAN=CH2_MINUS,ǰ���ͱ����ĻҶȾ�ֵ��
%ͬʱ�����������ڸ����������ǰ�����ظ���CH2_AREA�ͱ������ظ���CH2_BKD_AREA

%CH2_MEAN_MINUS=zeros(hnum-1,vnum-1);%ǰ���ͱ����ĻҶȾ�ֵ��
CH2_RAW=zeros(hnum-1,vnum-1);%ǰ���������صĻҶ�ֵ��
CH2_BKD_RAW=zeros(hnum-1,vnum-1);%�����������صĻҶ�ֵ��
CH2_AREA=zeros(hnum-1,vnum-1);%��¼(i,j)������ ǰ�����ص����
CH2_BKD_AREA=zeros(hnum-1,vnum-1);%��¼(i,j)������ �������ص����
CH2_MEAN=zeros(hnum-1,vnum-1);%ǰ���ĻҶȾ�ֵ
CH2_BKD_MEAN=zeros(hnum-1,vnum-1);%�����ĻҶȾ�ֵ

for i=1:hnum-1
    for j=1:vnum-1
        for ii=para2(i):para2(i+1)
         for jj=para3(j):para3(j+1)
%        for ii=para4(2*i)-t+1:para4(2*i+1)+t-1
%           for jj=para5(2*j)-t+1:para5(2*j+1)+t-1
%if ii<=h&&jj<=w
            if (III22(ii,jj)==1)
               CH2_RAW(i,j)=CH2_RAW(i,j)+PP2(ii,jj);%ǰ���ĻҶ�ֵ��
               CH2_AREA(i,j)=CH2_AREA(i,j)+1;%ǰ�������ظ���
            else
            
            CH2_BKD_RAW(i,j)=CH2_BKD_RAW(i,j)+PP2(ii,jj);%�����ĻҶ�ֵ��
            CH2_BKD_AREA(i,j)=CH2_BKD_AREA(i,j)+1;%���������ظ���
            end
%end
          end
       end
if CH2_BKD_AREA(i,j)==0    
   CH2_BKD_MEAN(i,j)=0;%�󱳾��ĻҶȾ�ֵ
else
   CH2_BKD_MEAN(i,j)=CH2_BKD_RAW(i,j)/CH2_BKD_AREA(i,j);%�󱳾��ĻҶȾ�ֵ  
end
% CH2_BKD_MEAN(i,j)=CH1_BKD_RAW(i,j)/CH2_BKD_AREA(i,j);%�󱳾��ĻҶȾ�ֵ
if CH2_AREA(i,j)==0  
   CH2_RAW(i,j)=1; 
   CH2_MEAN(i,j)=0;%��ǰ���ĻҶȾ�ֵ
        
else
   CH2_MEAN(i,j)=CH2_RAW(i,j)/CH2_AREA(i,j);%��ǰ���ĻҶȾ�ֵ
  
end
CH2_MEAN(i,j)=CH2_MEAN(i,j)-CH2_BKD_MEAN(i,j);
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%��ɹ���
flag=0;
medi2=median(median(CH2_AREA));
for i=1:hnum-1
    for j=1:vnum-1
        if CH2_AREA(i,j)>medi2*0.4&&CH2_AREA(i,j)<medi2*1.5
            flag=flag+1;
        end
    end
end
success2=flag/(i*j)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ʹ�þ���Ĵ�С��ͬ
% [H1,W1]=size(CH1_RAW);
% [H2,W2]=size(CH2_RAW);
% if H1>=H2
%     CH1_RAW=CH1_RAW(1:H2,1:W1);
% else  CH2_RAW=CH2_RAW(1:H1,1:W2);
% end
% [H1,W1]=size(CH1_RAW);
% [H2,W2]=size(CH2_RAW);
% if W1>=W2
%     CH1_RAW=CH1_RAW(1:H1,1:W2);
% else  CH2_RAW=CH2_RAW(1:H2,1:W1);
% end

%��log2ͼ��������ֵͼ
kmeans_VALUE=zeros((hnum-1)*(vnum-1));
% ch1=max(max(CH1_AREA));
% ch2=max(max(CH2_AREA));
for i=1:hnum-1
    for j=1:vnum-1
        if CH1_AREA(i,j)<10||CH2_AREA(i,j)<10||abs(log2(CH2_MEAN(i,j)/CH1_MEAN(i,j)))>10
           kmeans_VALUE(i,j)=10;
        else
           kmeans_VALUE(i,j)=log2(CH2_MEAN(i,j)/CH1_MEAN(i,j));
        end
    end
end
[HH,WW]=size(kmeans_VALUE);
%�������ֵ��0�ķ����ͣ�������SVM_VAR����
%�������ֵ��0��ƫ���ͣ�������SVM_DEV����
sum1=0;
sum2=0;
for i=1:HH
    for j=1:WW
        X=[kmeans_VALUE(i,j),0];
        sum1=sum1+(std(X))^2;
        sum2=sum2+abs(kmeans_VALUE(i,j));
    end
end

kmeans_VAR=sum1
kmeans_DEV=sum2

figure;
ii=1;
% hold on 
% axis([0 270 -8 10])
for i=1:(hnum-1)
    for j=1:(vnum-1)
       hold on    
        plot(ii,kmeans_VALUE(i,j),'b*');
%         title('ֱ��ͼ��ֵ�����ñ���ֵ');
        ii=ii+1;
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%��������ȡ��ı������ݶ�תΪexcel��ʽ������Ա�
xlswrite('k_CH1_AREA.xls',CH1_AREA)
xlswrite('k_CH1_BKD_AREA.xls',CH1_BKD_AREA)
xlswrite('k_CH1_MEAN.xls',CH1_MEAN)
xlswrite('k_CH1_BKD_MEAN.xls',CH1_BKD_MEAN)
xlswrite('k_CH1_RAW.xls',CH1_RAW)
xlswrite('k_CH1_BKD_RAW.xls',CH1_BKD_RAW)

xlswrite('k_CH2_AREA.xls',CH2_AREA)
xlswrite('k_CH2_BKD_AREA.xls',CH2_BKD_AREA)
xlswrite('k_CH2_MEAN.xls',CH2_MEAN)
xlswrite('k_CH2_BKD_MEAN.xls',CH2_BKD_MEAN)
xlswrite('k_CH2_RAW.xls',CH2_RAW)
xlswrite('k_CH2_BKD_RAW.xls',CH2_BKD_RAW)


xlswrite('k_VALUE.xls',kmeans_VALUE)
% xlswrite('k_VAR.xls',kmeans_VAR)
% xlswrite('k_DEV.xls',kmeans_DEV)
% t_k5=cputime-t_k5
% success1
% success2
end
